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Abstract 

An important part of building mathematical models 
based on measured data is calculating the accuracy 
associated with statistical estimates of the model 
parameters. Indeed, without some idea of this accuracy, 
the parameter estimates themselves have limited value. In 
this work, an expression for computing quantitatively 
correct parameter accuracy measures for maximum 
likelihood parameter estimates with colored residuals is 
developed and validated. This result is important because 
experience in analyzing flight test data reveals that the 
output residuals from maximum likelihood estimation are 
almost always colored. The calculations involved can be 
appended to conventional maximum likelihood estimation 
algorithms. Monte Carlo simulation runs were used to 
show that parameter accuracy measures from the new 
technique accurately reflect the quality of the parameter 
estimates from maximum likelihood estimation without the 
need for correction factors or frequency domain analysis of 
the output residuals. The technique was applied to flight 
test data from repeated maneuvers flown on the F-18 High 
Alpha Research Vehicle (HARV). As in the simulated 
cases, parameter accuracy measures from the new technique 
were in agreement with the scatter in the parameter 
estimates from repeated maneuvers, while conventional 
parameter accuracy measures were optimistic. 

Nomenclature 

a z vertical acceleration, g units 

c mean aerodynamic chord, ft 

C L lift coefficient 

D dispersion matrix 

E{ ■ } expected value 
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g gravitational acceleration, 32. 1 74 ft/sec 1 2 

l yy pitch axis moment of inertia, slug-ft 2 
J cost function 

m mass, slugs 

M information matrix 

n 0 number of outputs 

n p number of parameters 

N total number of sample times 

q body axis pitch rate, rad/sec 

q dynamic pressure, lbf/ft 2 

R discrete noise covariance matrix 

5R** autocorrelation matrix of vector v 

s sample standard error 

5 wing area, ft 2 

S(0 output sensitivity matrix at time (i - l)Ar 
t time, sec 

u(/) control vector 

V airspeed, ft/sec 

v(0 output residual vector at time (t - l)Af 
\(t) state vector 

y (t) output vector 

z(0 measured output at time (f - l)Ar 
a angle of attack, rad 

Sq Kronecker delta 

5 S stabilator deflection, rad 

At sample time, sec 

O roll angle, rad 

0 pitch angle, rad 

0 parameter vector 

V 0 gradient with respect to 0 

6 element of the parameter vector 

a Cramer-Rao bound for the standard error of 6 

\)(0 noise vector at time (/ - l)Ar 
0 zero vector 

V for every 
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superscripts 
T transpose 

estimate 
mean value 
time derivative 
-1 matrix inverse 

subscripts 

m measured 

o initial 

c corrected 

w wind axes 

Introduction 

Aircraft dynamic models typically include parameters 
which quantify the dependence of aerodynamic forces and 
moments on state and control variables. The values of 
these parameters are often estimated from flight test data. A 
good quantitative assessment of the accuracy of these 
parameter estimates is important for a variety of reasons 
related to experiment design, modeling, simulation, and 
flight control. 

Maximum likelihood 1 is commonly used to estimate 
aerodynamic parameters from flight test data. Assuming 
the model structure is correct, maximum likelihood 
parameter estimates approach the true parameter values, and 
the parameter variances approach their theoretical minimum 
values (the Cramer-Rao lower bounds), as the number of 
measured data points increases. Generally, a flight test data 
record length at least 2-3 times the period of the slowest 
dynamic mode to be modeled is sufficient for the parameter 
biases to be small and for the parameter variances to closely 
approach the Cramer-Rao bounds 2 * . In such cases, the 
Cramer-Rao bound can be used as a good approximation to 
the variance of maximum likelihood parameter estimates. 
References [2]-[4] compare and contrast the Cramer-Rao 
bound with other methods for assessing the accuracy of 
parameter estimates. Theoretical properties of maximum 
likelihood estimators and related arguments discussed in 
reference [2] indicate that the Cramer-Rao bound is the best 
accuracy measure for maximum likelihood parameter 
estimates. 

The research described here focuses on the output error 
formulation of maximum likelihood parameter estimation. 
This formulation includes measurement noise, but no 
process noise 1 . A modified Newton-Raphson optimization 
procedure 1 * 5 was used to determine the maximum 
likelihood parameter estimates. With this approach, the 
Cramer-Rao bounds are computed as part of the estimation 
procedure. It is well known, however, that the 
Cram6r-Rao bounds computed in this way are usually 
optimistic (too small) compared to the scatter in the 
parameter estimates from repeated flight test maneuvers 2 * 6 . 
This prompted the work of Maine, Iliff and 


Balakrishnan 1 * 2 * 7 - 9 , who traced the discrepancy to the fact 
that the residuals are colored for real flight test data analysis 
because some deterministic modeling error is always 
present. Output error techniques lump the deterministic 
modeling error together with the broad band random part of 
a measured signal and call this the measurement noise. 
This means the measurement noise is model -dependent and 
colored, because the modeling error usually lies in the same 
frequency band as the aircraft rigid body dynamics and 
accounts for a large part of the total noise power. 
References [l],[2],[7]-[9] describe how this kind of 
colored measurement noise is responsible for the 
discrepancy between the conventional calculation of the 
Cramer-Rao bounds and the observed scatter in 
flight-determined parameter estimates from repeated 
maneuvers. 

The theory underlying the output error formulation of 
maximum likelihood estimation assumes that the 
measurement noise is white Gaussian and band limited by 
the Nyquist frequency. The band limit is the result of 
discrete measurements taken at the sampling frequency, 
which is twice the Nyquist frequency. This measurement 
noise is broad band and incoherent. The term incoherent 
implies amplitude discontinuity and a lack of consistent 
phase-amplitude relationships, causing the autocorrelation 
function to be close to the impulse function. This part of 
the residual would be commonly recognized as having no 
deterministic component. If the structure of the model were 
correct, the residuals would be expected to be reasonably 
close to this type of noise. In real flight test data analysis, 
however, the residuals contain deterministic components 
from such sources as approximations to real aircraft 
aerodynamic dependencies, unmodeled dynamics such as 
structural modes, and linearization of the nonlinear 
equations of motion. The result is colored residuals which 
violate the assumptions of conventional maximum 
likelihood theory and lead to the aforementioned 
discrepancy 1 * 2 * 7 ' 9 . 

In reference [2], several engineering solutions were 
proposed to correct for the discrepancy. Each solution was 
based on the assumption that most of the residual power for 
real flight data analysis is concentrated in roughly the same 
frequency band as the rigid body dynamics and is due to 
deterministic modeling error. This assumption is stretched 
when relatively high frequency structural modes appear in 
the data or when the broad band random noise has a large 
enough magnitude to rival the power of the narrow band 
noise due to modeling error. For multiple outputs, the 
noise power from broad band random noise compared to 
that from narrow band deterministic modeling error is 
different for each output because of differences in the 
sensor characteristics and the physical quantity being 
measured. The solutions offered in reference [2] depend on 
knowing something about the bandwidth of the dominant 
source of power in the residuals. Obtaining this 
information requires Fourier transforms of the residuals and 
analysis in the frequency domain. The spectra of the 
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residuals depend on the model structure, the maneuver, the 
flight condition, and the instrumentation characteristics. All 
of these factors can change over the course of a flight test 
program, requiring changes in the corrections for the 
Cramer-Rao bounds. In addition, the band widths of the 
deterministic modeling error for the various measured 
outputs can be different from one another for the same 
maneuver. The solutions from reference [2] require some 
engineering judgment in the form of determining a 
correction factor or estimating the bandwidth of the 
dominant power in the residuals. Both of these approaches 
require an experienced analyst and limit the quantitative 
accuracy of the results. 

In the present work, a technique first put forth in 
reference [10] was used to process the residuals of a 
conventional maximum likelihood estimation in order to 
compute quantitatively accurate Cramer-Rao lower bounds 
for colored residuals. The approach accounts for colored 
residuals using a simple estimate of the residual correlation 
in the time domain. Existing maximum likelihood 
estimation routines can be easily upgraded because the 
technique involves a post-processing of the output residuals 
to correct the Cramdr-Rao lower bounds from the 
conventional calculation. The purpose of the present work 
is to document a few refinements to the technique described 
in reference [10], to validate the technique using Monte 
Carlo simulation with colored measurement noise, and to 
apply the technique to real data from repeated flight test 
maneuvers. 

The next section contains the theoretical analysis. 
Following this, the technique was applied in a controlled 
situation using simulated data from a model of the 
longitudinal dynamics of a fighter aircraft. The true 
parameter values were known, and the measured outputs 
were corrupted with colored noise, including both narrow 
band modeling error and broad band random noise. Using 
200 Monte Carlo simulation runs with various colored 
noise characteristics, it was demonstrated that the new 
technique produces Cramer-Rao bounds representative of 
the observed scatter in the estimates. The conventional 
Cramer-Rao bounds were found to be optimistic, in 
agreement with the results of previous research 1 ’2,6-9 

Next, the technique was applied to repeated longitudinal 
flight test maneuvers at 20 deg angle of attack for the F-18 
High Alpha Research Vehicle (HARV) aircraft. The scatter 
in the model parameter estimates from this flight test data 
was consistent with the Cramer-Rao bounds computed 
using the new technique, while the conventional calculation 
again gave optimistic values for the Cramer-Rao bounds. 


Theoretical Development 

The aircraft dynamic model can be represented as 

x(t) = f{x(t),u(t),B) (1) 

x(0j = x„ (2) 

y(t) = g(x(t), u(/j,0) (3) 

z(i) = y(i) + M<i) * = 1.2, N (4) 

The notation y(i) represents the sampled value of y(/) at 
f = (i-l)Af. There are N sampled data points. For 
conventional maximum likelihood, the discrete 
measurement noise vector V(i) is assumed to be zero mean 
white Gaussian and band limited at the Nyquist frequency, 

£{d<T)} = 0 £{u(j)« r O'j} = R<5, 7 (5) 

The maximum likelihood estimate of the parameter 
vector maximizes the conditional probability of realizing the 
measurements 1 ’ 5 : 

0 = arg max [ P( Z | 0 )] (6) 

where Z is the set of all measurement vectors z (/), for 
i=l, 2, N. The conditional probability distribution, 
P(Z |0) , also known as the likelihood function, is given 
by 


/> (z|e) = [(2^)"°|R|]" 


tex p\- 2 ^l z(i) ~y (i n R 1 [ z (0-y(i)] 


Maximizing the likelihood function in Eq. (7) is equivalent 
to minimizing its negative logarithm, known as the log 
likelihood function, 


-ln[p(Z|0)]= ~^[z(i)-y(i)] T R 1 [z (i)-y(i)\ 


+ y ln|R| 
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where the added constant 


n 0 N 

2 


\t\(2k) was omitted because 


it has no effect on the optimization. When R is known, 
minimizing the log likelihood function in Eq. (8) is 
equivalent to minimizing the cost function 


1 N 

J ( Q ) = - X [*<*> - y( v] T Rl [ z <‘) - y(*>] 


2 U 


(9) 


The cost in Eq. (9) can be minimized using a modified 
Newton-Raphson technique to determine parameter 
updates 1 ’ 4 5 , starting from some initial guess of the parameter 
vector. The initial guess for the parameter vector can be 
obtained from equation error methods 6 . Define the 
sensitivity matrix 


S (i) = 


dy(i)\ 


ae 


le=e 


* = 1 , 2 ,..., A (10) 


where theyth column of the sensitivity matrix contains the 
output sensitivities for the yth parameter, computed from 
central finite differences in Eqs. (l)-(3). The modified 
Newton-Raphson parameter update is given by 1510 


Often only the diagonal elements of the R matrix are 
estimated from Eq. (12), enforcing an assumption that the 
measurement noise sequences for the measured outputs are 
uncorrelated with one another. This assumption is 
generally a good one for real flight test data. All estimates 
of the measurement noise covariance matrix in this work 
assume a diagonal R matrix. Retaining the full R matrix 
could have been done with little conceptual difficulty, but 
the expected gains did not warrant the extra computation 
involved. The noise covariance matrix estimate, R, was 
used in the cost function of Eq. (9), and the minimization 
process described above for known R was repeated. Thus, 
the maximum likelihood estimation proceeds by alternately 
estimating the noise covariance matrix from Eq. (12) and 
minimizing the cost function using Eq. (11) with the latest 
value of the estimated noise covariance matrix. 
Convergence is reached when the estimated parameter 
vector 0, the estimated noise covariance matrix R , and the 
cost reach nearly constant values. Since maximum 
likelihood estimation is asymptotically unbiased 12 , the 
estimated parameter vector 0 should be close to the true 
value 0, and the gradient of the cost function with respect 
to the parameter vector should be close to zero. From Eq. 
(9), assuming R is held fixed. 


A0 = 0 — 0 = 


ZsW r R- 


S(i) 


i=l 


-1 


XSfi^R 1 [zft)-y(i)] 


v e^(0)|e = e = "I S(i) T R~ l [x(i)-y(i)) (13) 


1 = 1 


(ID 

The parameter vector update from Eq. (1 1) is added to 
the current estimate of the parameter vector in order to 
approach the true value of the parameter vector. In practice, 
there are times when the parameter vector update computed 
from Eq. (1 1) leads to an increase in the cost function or a 
divergence. This is because the modified Newton-Raphson 
step assumes that the current estimate of the parameter 
vector is near the true value. Using several iterations of a 
simplex algorithm 11 when the modified Newton-Raphson 
step produced an increase in the cost was found to be very 
effective in avoiding divergence and reaching a solution. 
This approach was followed in the present study. 

When repeated application of Eq. (11) converges, an 
estimate of the measurement noise covariance matrix, R, 
can be obtained from the output residuals. The expression 
for the estimate of R comes from taking the derivative of 
the right hand side of Eq. (8) with respect to R, setting the 
result equal to zero, and solving for R, 

1 N 

( 12 ) 


For practical computation, simultaneous satisfaction of 
the numerical criteria given below were used to define 
convergence of the maximum likelihood estimation: 


im-IH-, 

fc-L -faLi 


IH-i 



i-yf 

0*-.) 

J 

(®*-i 

) 


<1.0x10-5 V j, ; = 1,2 n. 


<0.05 Vi, 1 = 1,2,... 


< 0.001 


(14) 


a/(e) 


ae. 


<0.05 V j, j = l,2,...,n„ 


J b=e 


where k denotes the current estimate iteration number and 
n, denotes the estimate of the ith diagonal element of R. 
The approximate expression for the cost gradient with 
respect to the parameters (Eq. (13)) was used for the last 
criterion in (14). 

The minimum achievable parameter variances, the 
Cramer-Rao lower bounds, are given by the diagonal 
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elements of the dispersion matrix, D ^2,5 This dispersion 
matrix is defined as the inverse of the information matrix 
M, the latter being a measure of the information contained 
in the data from an experiment. The expressions for these 
matrices are 1 - 2 5 


= D 


XX SO) r R-^v( ( >or}R-‘SO-J 


.1=1 7*1 


D 


(20) 


N 

M = ^S(i) T R^S(i) (15) 

i=l 


D = M -1 = 


N 

I 

L>i 


-i-i 


S(i) T R~' S (i) 


(16) 


The square root of the yth diagonal element of D, djj, 
gives the Cram6r-Rao lower bound for the standard error of 
the /th parameter estimate, 


a i=^ij j = x > 2 n P ( 17 > 


It can be seen from Eqs. (11) and (16) that the 
dispersion matrix is computed when determining the 
modified Newton-Raphson step as part of the conventional 
maximum likelihood estimation . The assumption that the 
output residuals are white and therefore uncorrelated in time 
is implicit in the algorithm and indicated in Eq. (5). The 
next section details the theory involved in accounting for 
arbitrary colored output residuals, which are correlated in 
time. 

When the conventional maximum likelihood estimation 
has converged, the estimated parameter vector will be close 
to the true value and Eq. (11) holds. Define the residual 
vector 


\(i) = z(i)-y(i) i = l, (18) 


When the output residuals are assumed to be zero mean 
white (cf Eq. (5)), then 


E{\(i)v(j) T } = RS iJ (21) 


From Eqs. (16), (20) and (21) it is easy to see that the 
parameter covariance matrix reduces to the dispersion 
matrix D when the output residuals are white. 

For colored residuals, 


£{v(/)vO) r } = SR vv (/-y) (22) 


where is the autocorrelation of the output 

residuals, so that the estimated parameter covariance matrix 
can be computed from Eq. (20) using an estimated value for 
9*w0' ~ 7 ). Define the discrete unbiased estimate of the 
output residual autocorrelation 12 

1 N ~ k 

*w(*) = ^v(i)v(i + k) T = (23) 

1 = 1 

Substituting Eq. (22) into Eq. (20) results in 


cov 


(e)-i> 


£S(/) r R-‘X *„(/-,•) R-'SO') 

.* = 1 7 = 1 


D 

(24) 


The estimated parameter covariance matrix can be expressed 
using Eq. (11) with substitutions from the definitions in 
Eqs. (16) and (18), 

cov(e) = £j(e-e)(e-e) r j 

f N N ) 

= £ XX DS ^ rR_1 y(i)y(j) T R- l S(j)D\ 


If it is assumed that the discrete noise covariance matrix 
and the output sensitivities do not depend on the parameter 
vector estimate at the maximum likelihood solution, then the 
estimated parameter covariance matrix can be written as 


Eq. (24) was used to account for colored residuals, 
which are correlated in time. Eq. (23) was used to estimate 
j) in Eq. (24). The values for D, R _1 , and S are 
from the conventional maximum likelihood estimation. 
Eqs. (23) and (24) embody the post-processing applied to a 
conventional maximum likelihood solution to account for 
colored residuals. 

If the colored output residuals are assumed to be caused 
by modeling error, then the maximum likelihood parameter 
estimates are biased, so that 

£{e} = e+b(e) (25) 

where b(0) is a vector of bias errors which are unknown 
and unknowable in practice. For these biased parameter 
estimates, the Cramer-Rao bound will contain terms in 
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addition to D, namely D multiplied by V e b(0) and 
[V e b(0)][V e b(0)] r (see reference [2] for details). In the 
present analysis, the deterministic modeling error is 
included as part of the measurement noise, making the 
residuals colored. The approach taken here for arbitrary 
frequency content in the measurement noise could then be 
viewed as effectively including the bias error due to 
deterministic modeling error as part of the parameter 
accuracy measure that accounts for colored residuals. 

Results 

The longitudinal short period dynamics of the F-18 
HARV fighter aircraft at approximately 20 degrees angle of 
attack were studied. The model state equations in wind 
axes are given by 


or 


= -^rrfcz a +C Z — + C Z 6. \ + q 
mV{ “ i 2V s * 1 1 

+ ~(cos( <P m ) cosf Q m ) cos( a) + sin( O m ) sinf a)) + a* 


with measurement equations 

a m (i) = a(i) + v x (i) 
<lm(‘) = q<i) + v 2 (i) 


a z Ji) 


= ( C z a(i) + C z + C z 8 (i) 

mg V ° i 2V 


(26) 


+ a, + v 3 (i) 


i = 1,2, .... N 
(27) 


assuming that C Z =-C L , a lv = a z , and where 

«„* = (qS/mV) C Zo + d 0 and q* = (qSc/l„) C Mg + 4o ■ 

Initial conditions for the states were computed from the 
measured time histories of a and q using a time domain 
smoother 13 . 

To validate the new technique for computing 
Cram6r-Rao bounds, two hundred Monte Carlo simulation 
runs were made using various colored measurement noise 
processes. Each noise sequence had part of its power in the 
frequencies between 0 hz and 1 hz inclusive (roughly the 
frequency band of the uncorrupted simulation outputs), 
with the remaining power taken up by white Gaussian noise 
out to the Nyquist frequency. The narrow band portion of 
the colored noise sequence was generated by passing zero 
mean white Gaussian noise through a fifth order 


Chebyshev low pass filter with frequency cutoff set at 
1 hz. The resulting narrow band noise was combined with 
wide band noise from a separate realization of the zero 
mean white Gaussian noise process. The percentage of the 
total noise power from the narrow band noise was 
determined for each Monte Carlo run by a random number 
with uniform distribution on the interval [0,100]. This 
procedure was carried out for each simulated output on each 
Monte Carlo run, and the resulting colored noise sequences 
were scaled to achieve approximately a 5 to 1 signal to 
noise ratio for all simulated outputs. Figure 1 shows the 
power spectral density for the colored noise added to a for 
run 200, where 19% of the noise power was in the 
frequency range of 0 hz to 1 hz, inclusive. Colored noise 
sequences generated in this way are representative of 
residual sequences observed when analyzing real flight test 
data, and were used for that reason. 

To make the Monte Carlo simulation runs realistic, the 
stabilator input was taken from measured data for the F-18 
HARV flying a maneuver designed specifically for accurate 
parameter estimation 14 . The stabilator input is shown as 
the solid line in figure 2. The values of the parameters used 
in the simulations (given in column 2 of Table 1) 
approximately reflect the short period dynamics of the F-18 
HARV at 20 degrees angle of attack. The stabilator input 
and parameter values were the same for each simulated data 
run, so that the information in the data was constant. The 
sampling rate was 50 hz and the data record length was 14 
seconds. Maximum likelihood estimation as described in 
the previous section was used to estimate the parameters. 

Since the parameter values were known for the 
simulated data, the true accuracy of the maximum likelihood 
estimates could be compared to the accuracy indicated by 
the Cramer-Rao bound calculations. The conventional 
Cramer-Rao bounds for the parameter standard errors were 
denoted by <7 and were computed as the square root of the 
diagonal elements of matrix D in Eq. (16). The 
Cramer-Rao bounds for the parameter standard errors 
corrected for colored residuals were denoted by cr c and 
were computed as the square root of the diagonal elements 
of the covariance matrix from Eq. (24). Results from both 
the conventional computation and the corrected calculation 
were expressed in terms of the ratio of the absolute 
deviation of each parameter estimate from its true value to 
the computed Cramer-Rao bound for the parameter standard 
error. This accuracy measure was assigned the symbol Tf y 


r? = [ 0 - 0 1 /<r 
= | 0 — 0 1 /o c 


(28) 


where the subscript c denotes values for the corrected 
Cram6r-Rao bounds. 

For a maximum likelihood estimator, the probability 
distribution of the parameter estimates about their true value 
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approaches a Gaussian distribution as the number of data 
points gets large. Evidence of this can be found in figure 3, 
which is a histogram of the parameter estimates from all 
two hundred Monte Carlo runs for the C M a parameter. 

Corresponding histograms for the other estimated 
parameters were similar. It follows that 7] and rj c should be 
less than 3 almost all the time if the computed Cramer-Rao 
bounds reflect the true accuracy of the parameter estimates. 

Table 1 shows results for two representative Monte 
Carlo runs. Columns 4 and 5 for run 47 and columns 7 
and 8 for run 185 show that the corrected Cramer-Rao 
bounds accurately reflected the true parameter accuracy, 
while the conventional Cramer-Rao bounds were optimistic 
(i.e., too small) and produced 7] ratios which exceeded 3 
for almost every estimated parameter. Considering the full 
set of two hundred Monte Carlo runs, Table 2 gives the 
mean values and standard errors of 7] and rj c for each 
estimated parameter. This data shows that the conventional 
Cramer-Rao bounds were inaccurate on the average and 
exhibited a large variability, while the converse was true for 
the corrected Cramer-Rao bounds. 

Table 3 gives another summary of the parameter 
estimation results for the two hundred Monte Carlo 
simulation runs. The second column of the table gives the 
mean values of the parameter estimates, and the third 
column gives the sample standard errors for the parameter 
estimates, s, computed from the scatter of the parameter 
estimates. Columns 4 and 6 give the mean values of the 
Cramer-Rao bounds for the parameter standard errors 
computed using the conventional and corrected 
techniques, a and a c , respectively. Columns 5 and 7 
show the ratio of the sample standard errors for the 
parameter estimates to a and a c , respectively. These 
values are far less than 3 for the corrected calculation of the 
Cramer-Rao bounds, indicating a proper accounting for the 
changes in the residual spectra, while the conventional 
calculation of the Cramer-Rao bound was optimistic, 
producing values of the s/a ratio greater than 3. 

The data in Tables 1,2, and 3 demonstrate that the 
extent to which the conventional Cramer-Rao bounds 
misrepresented the true parameter accuracy was neither 
consistent nor predictable from parameter to parameter or 
from run to run. This phenomena has been observed 
previously when analyzing flight test data from repeated 
maneuvers 6 . It follows that the common practice of 
applying a fixed correction factor to the conventional 
calculation of the Cramer-Rao bounds is incorrect to a 
varying and unpredictable degree in cases where coloring of 
the residual spectrum varies, as in this simulation study. 
Changes in the coloring of the residuals similar to those 
studied here can easily be brought about in practice by 
changes in the model structure, the maneuver, the flight 
condition, or the instrumentation. 

Next, flight test data was analyzed from five repeats of 
the same longitudinal maneuver, flown on the F-18 HARV 
at approximately 20 degrees angle of attack and 25,000 


feet. The input was applied to the symmetric stabilator by a 
computerized On-Board Excitation System (OBES), so that 
the runs were very nearly repeats of one another. Figure 2 
shows the excellent repeatability using the OBES system 
for the five repeated flight test runs of the stabilator input 
maneuver. All of the data used for analysis was sampled at 
50 hz. Corrections were applied to the angle of attack and 
accelerometer measurements to account for sensor offsets 
from the center of gravity, and the angle of attack 
measurement was corrected for up wash. Data compatibility 
analysis 15 revealed that the data from the sensors was 
consistent to a degree which warranted no further 
corrections. Maximum likelihood estimation was carried 
out using the same procedure as for the Monte Carlo 
simulation runs. The same model given in Eqs. (26)-(27) 
was used. 

Table 4 gives flight test results in a form similar to 
Table 3. Column 7 shows that the corrected Cramer-Rao 
bounds were an accurate measure of the scatter in the 
parameter estimates. In column 5, the conventional 
Cramer-Rao bounds were again optimistic for the pitching 
moment (C M ) parameters but were close to correct for the 
vertical force ( C z ) parameters. The reason is that the a and 
a z measurements are the main influences on the C z 
parameters, and the residuals for both these outputs 
exhibited considerable power at high frequencies, due to 
unmodeled structural modes. The power spectrum for a 
typical a z residual (from run 1) is shown in figure 4. 
These colored residual spectra roughly resembled constant 
power out to the Nyquist frequency, which is the 
assumption made in the theory underlying the conventional 
Cramdr-Rao bound calculation. The q measurement did not 
have these high frequency components, so the conventional 
Cram6r-Rao bound calculation gave very optimistic values 
for the Q/ parameters. The power spectrum for a typical q 
residual (from run 1) is shown in figure 5. The corrected 
calculation of the Cram6r-Rao bound worked equally well 
for the C z and Cm parameters because information about 
the particular coloring of the residuals is incorporated 
automatically via the autocorrelation estimate from Eq. (23) 
used in Eq. (24). 

Figures 6 and 7 depict the parameter estimation results 
for the aerodynamic parameters. The error bars represent 
the Cramer-Rao bounds for the standard errors computed 
using the conventional calculation for figure 6 and the 
corrected calculation for figure 7. The error bars in the 
lower three plots of figure 6 are difficult to see because they 
are roughly the size of the circles on the plot marking the 
individual parameter estimates. These plots and the 
accompanying data in Table 4 show that the standard 
calculation for the Cramer-Rao bounds gave optimistic 
values compared to the scatter in the estimates from 
repeated maneuvers, whereas the corrected calculation for 
the Cramdr-Rao bounds produced Cramer-Rao bounds 
which accurately reflected the scatter of the estimates. 


7 

American Institute of Aeronautics and Astronautics 



Concluding Remarks 

Algorithms for aircraft parameter estimation using the 
output error formulation of maximum likelihood are in 
widespread use. The Cramer-Rao bounds characterizing 
parameter accuracy that are obtained from conventional 
calculations are known to be generally optimistic (i.e., too 
small) in practice, compared to the scatter in parameter 
values estimated from repeated maneuvers. Parameter 
estimates have limited utility when there is no firm idea of 
their accuracy. In this work, an expression for correcting 
Cram£r-Rao bounds from maximum likelihood estimation 
with colored residuals was presented and validated. This 
result is important because the residuals from maximum 
likelihood estimation are almost always colored in practice, 
due to deterministic modeling error. 

The calculations involved in the algorithm for 
computing Cramer-Rao bounds that account for colored 
residuals can be carried out in a short subroutine called at 
the conclusion of a conventional maximum likelihood 
estimation algorithm. Bandwidth of the dominant power in 
the residuals need not be known or estimated, since it is 
accounted for automatically in the algorithm by an unbiased 
estimate of the residual autocorrelation. There is no need 
for correction factors. The algorithm was shown to work 
for a wide range of colored residual spectra similar to what 
might be encountered in real flight test data analysis. All 
calculations are performed in the time domain, obviating the 
need for frequency domain analysis of the residuals. 

The corrected calculation for the Cramdr-Rao bounds 
presented here produced consistently accurate measures of 
the scatter in the parameter estimates, using an algorithm 
with moderate computational cost that can be applied as a 
post-processing of the output residuals from a conventional 
maximum likelihood solution. 

Monte Carlo simulation runs using various colored 
noise sequences were carried out to validate the 
performance of the algorithm. Analysis of flight data from 
repeated maneuvers flown on the F-18 High Alpha 
Research Vehicle (HARV) demonstrated the validity of the 
technique for computing appropriate parameter accuracy 
measures using real flight test data. The algorithm 
described in this work was shown to be an effective means 
of accurately determining the quality of parameter estimates 
from the output error formulation of maximum likelihood 
estimation. 
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Table 1 PARAMETER ESTIMATION RESULTS Table 3 PARAMETER ESTIMATION RESULTS 

FOR TYPICAL MONTE CARLO RUNS FOR 200 MONTE CARLO RUNS 



e 


Run 47 


Run 185 


Simulation 

Conventional 

Corrected 

0 

n 

V c 

6 

n 

4c 

0 

s 

a 

s/<7 


s/a c 

c z 

z. a 

>2.00 

-1.30 

4.24 

0.94 

-3.78 

9.33 

2.26 

C z 

t'CC 

-1.94 

0.819 

0.179 

4.56 

0.521 

1.57 



-13.8 

4.10 

0.99 

-81.8 

1.37 

0.33 

c z. 

-62.1 

51.1 

12.3 

4.17 

38.7 

1.32 

% 

-0.90 

-1.48 

4.20 

0.95 

0.04 

6.46 

1.59 

Cz * s 

-0.92 

0.550 

0.139 

3.96 

0.441 

1.25 

“o 

-0.80 

-1.08 

4.73 

1.04 

0.01 

10.46 

2.33 


-0.82 

0.312 

0.068 

4.56 

0.199 

1.57 

C M a 

-0.30 

-0.35 

7.00 

1.67 

-0.32 

3.09 

0.48 

C M a 

-0.30 

0.037 

0.008 

4.80 

0.022 

1.67 

C M q 

-16.0 

-17.3 

3.55 

0.86 

-13.2 

7.04 

1.13 

C M, 

-16.1 

1.822 

0.448 

4.07 

1.278 

1.43 


-0.70 

-0.67 

5.42 

1.03 

-0.72 

3.86 

0.57 

CM Ss 

-0.70 

0.030 

0.007 

4.11 

0.021 

1.43 


0.08 

0.09 

5.77 

1.20 

0.09 

5.01 

0.73 

4*0 

0.08 

0.013 

0.003 

4.67 

0.008 

1.64 

a z 

z o 

-0.70 

-1.07 

6.11 

1.36 

0.04 

9.48 

2.14 

“to 

-0.72 

0.313 

0.069 

4.54 

0.198 

1.58 


Table 2 PARAMETER ACCURACY MEASURE 
STATISTICS FOR 200 MONTE CARLO RUNS 


Conventional Corrected 



n 


4c 

a n c 

Cz 

z a 

3.84 

2.85 

1.37 

1.10 

Cz, 

3.43 

2.72 

1.10 

0.83 

Cz Ss 

3.30 

2.51 

1.10 

0.91 

d* 

3.88 

2.89 

1.36 

1.05 

C M a 

4.03 

3.42 

1.44 

1.38 

C M g 

3.62 

2.66 

1.26 

1.08 

C M Ss 

3.63 

2.72 

1.23 

0.88 

4* 

3.97 

3.19 

1.37 

1.15 

a z 

z o 

3.88 

2.80 

1.39 

1.06 


Table 4 PARAMETER ESTIMATION RESULTS 
FROM FLIGHT TEST DATA 



Flight 

Conventional 

Corrected 


e 

s 

a 

s/a 

Oc 

s/a c . 

c z 

^ a 

-2.04 

0.161 

0.083 

1.95 

0.192 

0.84 

Cz « 

-61.0 

4.86 

3.48 

1.40 

7.67 

0.63 

Cz S s 

-0.92 

0.081 

0.044 

1.82 

0.081 

1.00 

«« 

-0.72 

0.068 

0.031 

2.21 

0.068 

1.00 

C M a 

-0.30 

0.065 

0.006 

10.80 

0.052 

1.26 

c m. 

-19.1 

2.55 

0.372 

6.86 

2.84 

0.90 

C m Ss 

-0.74 

0.048 

0.007 

6.90 

0.051 

0.94 

m 

0.08 

0.025 

0.002 

10.62 

0.020 

1.24 

cl. 

z o 

-0.66 

0.067 

0.031 

2.17 

0.070 

0.95 
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Figure 1 Example colored noise power spectrum, 
run 200 



Figure 4 Power spectrum of a z residuals, 
flight test run 1 



Figure 2 Stabilator input time histories from five 
repeated flight test maneuvers 
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Figure 5 Power spectrum of q residuals, 
flight test run 1 
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